package fem2.element;

import fem2.Edge;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.Node;
import fem2.NotImplementedError;
import fem2.enu.MassMatrixType;
import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IMatrixReadOnly;
import inf.jlinalg.IVector;
import inf.jlinalg.IVectorReadOnly;
import inf.math.LinAlg;
import math2.LinAlgUtilities;

public class TrussElement extends StructuralElement {

	private IMatrix T;

	public TrussElement(MeshPart mp, MaterialModel mm) {
		super(mp, mm);
		if (!(mp instanceof Edge)) {
			throw new Error("truss element must be defined on an edge");
		}
	}

	public Edge getEdge() {
		return (Edge) getMeshPart();
	}

	@Override
	public void computeB(MeshPart mp, IMatrix b, double... xi) {
		int n = mp.getBasis().getDimension();
		Edge e = getEdge();
		double[][] gGradients = e.gGradientsAt(xi);
		for (int i = 0; i < n; i++) {
			double g1 = gGradients[i][0];
			b.set(0, i, g1);
		}
	}

	@Override
	public IMatrixReadOnly computeKte(IMatrix ke) {
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		Array1DMatrix B = Array1DMatrix.createColumnMajor(1, n);
		double a = getMaterialModel().getSectionProperty();
		int ng = getNumIntegrationPoints();
		IMatrix kk = new Array2DMatrix(n, n);

		for (int k = 0; k < ng; k++) {
			double w = getIntegrationWeight(k);
			double[] xi = getIntegrationPoint(k);

			computeB(mp, B, xi);
			double[][] jac = mp.jacXAt(xi);
			double detJ = LinAlg.norm2(dim, jac[0]);
			double coeff = w * detJ * a;

			IMatrix C = getMaterialModel().getConstitutiveMatrix(this, k);
			LinAlgUtilities.btCB(B, C, coeff, kk);
		}

		LinAlgUtilities.btCB(T, kk, 1.0, ke);

		return ke;
	}

	@Override
	public IMatrixReadOnly computeKge(IMatrix ke) {
		/*
		 * TODO
		 */
		throw new NotImplementedError();
	}

	@Override
	public IMatrixReadOnly computeMe(IMatrix me, MassMatrixType massType) {
		/*
		 * TODO
		 */
		throw new NotImplementedError();
	}

	@Override
	public IVectorReadOnly computeRie(IVector rie) {
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		Array1DMatrix B = Array1DMatrix.createColumnMajor(1, n);
		double a = getMaterialModel().getSectionProperty();
		int ng = getNumIntegrationPoints();
		IVector rr = new ArrayVector(n);

		for (int k = 0; k < ng; k++) {
			double w = getIntegrationWeight(k);
			double[] xi = getIntegrationPoint(k);

			computeB(mp, B, xi);
			double[][] jac = mp.jacXAt(xi);
			double detJ = LinAlg.norm2(dim, jac[0]);
			double coeff = w * detJ * a;

			IVector s = getStress(k);
			BLAM.multiply(coeff, BLAM.TRANSPOSE, B, s, 1.0, rr);
		}

		BLAM.multiply(1.0, BLAM.TRANSPOSE, T, rr, 1.0, rie);

		return rie;
	}

	@Override
	public IVector computeRegularStrain(double... xi) {
		MeshPart mp = getMeshPart();
		IVector u = new ArrayVector(mp.getUHatLocal());
		int d = 1;
		int n = mp.getBasis().getDimension();
		int ntrs = 1;
		Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, d * n);
		computeB(mp, B, xi);
		IVector ulocal = new ArrayVector(n);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, T, u, 0.0, ulocal);
		IVector e = new ArrayVector(ntrs);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, B, ulocal, 0.0, e);
		return e;
	}

	@Override
	public IVectorReadOnly computeRes(IVector res) {
		/*
		 * TODO
		 */
		throw new NotImplementedError();
	}

	@Override
	public void Initialize() {
		/*
		 * initialize the sigma array for this element
		 */
		int ng = getNumIntegrationPoints();
		for (int i = 0; i < ng; i++) {
			sigma.add(new ArrayVector(1));
		}
		/*
		 * initialize material
		 */
		getMaterialModel().Initialize(this);
		/*
		 * compute rotation matrix
		 */
		MeshPart mp = getMeshPart();
		int dim = mp.getDimension();
		int n = mp.getBasis().getDimension();
		T = Array1DMatrix.createRowMajor(n, dim * n);
		Node[] nodes = mp.getNodes();
		double[] n1 = nodes[0].getPosition();
		double[] n2 = nodes[1].getPosition();
		double[] dd = new double[dim];
		for (int i = 0; i < dim; i++) {
			dd[i] = (n2[i] - n1[i]);
		}
		double L = LinAlg.norm2(dim, dd);
		double[] cos = new double[dim];
		for (int i = 0; i < dim; i++) {
			cos[i] = dd[i] / L;
		}
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < dim; j++) {
				T.set(i, dim * i + j, cos[j]);
			}
		}
	}

}
